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certain statistical distributions which are of importance in linguistics, 
statistical-mechanics theory, and biophysics (statistical analysis of 
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in experimental mathematics, and we briefly expand on further ap- 
plications in theoretical physics. Finally, we discuss a related Mathe- 
matica program for the computation of Lerch's transcendent. 
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1 Introduction 



The combined nonlinear-condensation transformation ||T| (CNCT) is an algorithm that trans- 
forms a sequence of partial sums of a slowly convergent nonalternating (monotone) series (i.e., 
a series whose terms all have the same sign) into a sequence of transforms (to be referred to as 
CNC transforms) often with better numerical properties. A number of applications (e.g. [§,^) 
have recently been described; an acceleration of the convergence by several orders of magni- 
tude is observed in many cases. The CNCT addresses, at least in part, the difficulties usually 
associated with the acceleration of the convergence of monotone series. The application of con- 
vergence acceleration methods to series of this type typically leads to numerical instabilities and 
cancellations in higher transformation orders (for a discussion see [^). In many cases, it has 
been observed that the CNCT is a remarkably stable numerical process, mainly because the 
nonalternating input series is transformed into an alternating series before the actual conver- 
gence acceleration method is applied. The CNCT is a process involving two etappes; first step: 
transformation nonalternating^ alternating series, second step: convergence acceleration via the 
delta transformation (this transformation was introduced in |5|, and its usefulness as a powerful 
generalized resummation and convergence acceleration method was demonstrated in Q). 
The purpose of this paper is threefold: (i) to discuss applications of the CNCT in such diverse 
disciplines as statistics, theoretical physics, and experimental mathematics, (ii) to propose an 
implementation of the CNCT for the calculation of Lerch's transcendent in Mathematica 
and (iii) to provide numerical evidence for the convergence of the CNCT in higher transformation 
order using computer arithmetic with enhanced accuracy |^-10|. 

Numerical evidence gained from high-precision calculations might be useful because, currently, 
no direct, general proof exists for the convergence of the delta transforms |Q, ||, which form 
part of the second step of the CNCT. Of course, it should be noted that the convergence of the 
delta transforms has been proven for a set of model problems j^, and the exactness of the CNC 
transformation has been shown in for a number of model series, but this information is of 
limited use for the investigation of realistic applications. This fact has to be contrasted with 
the observed rapid convergence of the transforms in many applications of practical importance 
(e.g. |§, ^). In particular, the recent evaluation Q of the bound-electron self-energy has lead 
to results with roughly 19 significant figures for atomic hydrogen (the CNCT has been used in 
this calculation for the acceleration of slowly convergent partial wave expansions). For atomic 
hydrogen, the renormalization process leads to a numerical loss of about 9 significant figures, 
because the physically significant part of the bound-electron self-energy is obtained after sub- 
tracting the self-energy of the free electron, the latter is a part of the electron mass. The succesful 
verification of the consistency of the remaining 10 significant figures left after the renormaliza- 
tion with higher-order analytic results |11] provides a sensitive test of the numerical results 
and, consequently, of the convergence properties and the reliability of the algorithm described 
in Q and in the current article. To complement this evidence, we investigate here the CNCT 
with extended-precision arithmetic and gauge the rate of convergence of the CNC transforms in 
higher order. 

This paper is organized as follows: in Sec. we recall the motivation and the formulas used 
in the construction of the CNCT. In Sec. ^, we discuss the application of the CNCT to the 
evaluation of a class of important special functions, which may be used e.g. for the statistical 
analysis of DNA sequences. In this context we also discuss an explicit implementation of the 
CNCT in Mathematica which may serve as a template for applications in other areas and pro- 
gramming languages. These investigations are supplemented by the high-precision calculations 
in Sec. ^ where a certain application is discussed which could probably best be associated to 
the field of "experimental mathematics" . In Sec. ^, we discuss further applications of the CNCT, 
mainly related to the evaluation of other special functions which occur naturally in theoretical 
physics. We provide concluding remarks in Sec. ^ Appendix ^ contains a detailed description 
of Mathematica code for the calculation of Lerch's transcendent. 



2 The Combined Nonlinear— Condensation Transformation 



2.1 CNCT: Motivation 

We briefly recall the notion of logarithmic convergence and the difficulties associated with the 
acceleration of the convergence of nonalternating series. Let a sequence {sn}'^=Q fulfil the asymp- 
totic condition 

hm = p, (1) 

n^oo Sn — S 

where s = Soo is the limit of the sequence as n — > oo. If p > 1, then the sequence {s^jJ^Q is 
divergent. For p = 1, the sequence may either be convergent or divergent. A convergent sequence 
with \p\ = 1 is called logarithmically convergent (if p < 1, the series is called linearly convergent). 
Let us further assume that the elements of the sequence {sn}^=o in Eq. (||) represent partial 
sums 



Y,a{k) (2) 



k=Q 

of an infinite series. Here, we will almost exclusively investigate slowly convergent nonalternating 
sequences {sn}'^=Q whose elements are all real and positive (for these sequences, < p < 1). In 
the case of slow convergence, p is either very close or equal to unity. 

As observed by many authors (e.g. [Q]), the acceleration of the convergence of nonalternating 
sequences is a potentially unstable numerical process. The reason is the following: A sequence 
transformation can only accelerate convergence if it succeeds in extracting additional information 
about the index-dependence of the remainders or "truncation errors" 

1"n — Sn ~ S (3) 

from a necessarily finite set of partial sums {sn}n=o input series. Normally, this is done by 

forming arithmetic expressions involving higher weighted differences of the s„. The calculation 
of higher weighted differences is a potentially unstable process which can easily lead to a serious 
loss of numerical significance due to cancellation if the input data all have the same sign. 

The main notion of convergence acceleration is to extract information "hidden in trailing digits" 
from a necessarily finite number of sequence elements, in order to convert a sequence {sn}^=Q 
into a new sequence {s^jj^g with hopefully better numerical properties. 

Let us assume that the input sequence {sn}'^=Q is logarithmically or linearly convergent (with 
p < 1) and that a sequence transformation can be constructed so that the transformed sequence 
{s^lJ^Q is linearly convergent with p' < p [see Eq. (|l]) for the definition of p\. Convergence of 
the new sequence can be said to be accelerated if p' < p. An equivalent definition of convergence 
acceleration can be given by calculating the limit of the following sequence of the ratios Xn, 

Xn = ^, (4) 
Sn S 

and to define convergence to be accelerated if 

lim Xn = 0- (5) 

n— *oo 

Indeed, if p' < p, then 

lim ^ EE hm ^^^±1^ ^li^ = ^ < 1. (6) 

Xn S'n- S Sn+l - S p 

So, if p' < p, then Xn asymptotically behaves as a geometric progression within its circle of 
convergence, and thus lim„^oo Xn = 0- 



2.2 CNCT: Formulas 



Recently, the combined nonlinear-condensation transformation [|| (CNCT) has been proposed 
as a computational tool for the accelerated numerical evaluation of slowly convergent, nonalter- 
nating series. The idea is to divide the acceleration process in two steps. The first step, which 
is a reordering process, consists in a rearrangement of the terms of the input series into an 



alternating series via a Van Wijngaarden transformation |12|. The output of the first step is an 
alternating series whose terms do not decay more rapidly in magnitude than those of the original 
input series (see the discussion in Sec. 3 of |]l|). It could appear that nothing substantial has 
been achieved in the first step of the CNCT. The Van Wijngaarden step of the CNCT merely 
represents a "computational investment" with the intention of transforming the nonalternating 
input series into a form which is more amenable to the acceleration of convergence. The second 
step, which represents a convergence acceleration process, consists in the application of a power- 
ful nonlinear sequence transformation for the acceleration of the convergence of the alternating 
series which resulted from the first step of the CNCT. 



Following Van Wijngaarden |12|, we transform the nonalternating input series 



oo 



J2a{k), a{k)>0, (7) 

whose partial sums are given by (H), into an alternating series '^'jLoi~^y -^j • After the first step 
of the transformation, the limit of the input series is recovered according to 

oo oo 

Y,a{k) = Y^{-iyAj. (8) 

k=0 j=0 

The quantities Aj are defined according to 

c 

A, = E , (9) 



k=0 



where 

bjf') = 2'= 0(2^= (i + 1) - 1) . (10) 

The Aj are referred to as the condensed series and the series Yl'jLoi~^y -^j referred to as 
the transformed alternating series, or alternatively as the Van Wijngaarden transformed series. 

The construction of the condensed series reminds one of Cauchy's condensation theorem (see 
e.g. p. 28 of Ref. or p. 121 of Ref. |l^). Given a nonalternating series YlT=o ^(^) with terms 
that satisfy \a{k + 1)| < |a(/!;)|, Cauchy's condensation theorem states that YlT=o "^(k) converges 
if and only if the first condensed series Aq defined according to Eq. (^) converges. 

The summation over k in Eq. (^) does not pose numerical problems. Specifically, it can be easily 
shown in many cases of practical importance that the convergence of Xlfe^o ^) linear 

even if the convergence of Yl^=o o,{k) is only logarithmic. We will illustrate this statement by 
way of two examples. Example 1: a logarithmically convergent input series whose terms behave 
asymptotically as a{k) ^ k~^~'^ with e > 0. In this case, the partial sums 

n 

Af ) = bi^-) (11) 



k=0 



converge linearly with 



A^"+^^-A- 1 
hm = , < 1 , a{k) ~ k-'-^ , - oo. (12) 

n^oo ^(") _ ji^. 2^ (j + 1)^+^ 



Example 2: a series with a{k) ~ k^r^ where < r < 1 and /? real. Here, we have p = r < 1, and 
the series is (formally) linearly convergent. However, slow convergence may result if p is close to 
one. In this case, the condensed series are very rapidly convergent, 

hm ^ = 0, a(A;)~A:V, k ^ oo. (13) 



n— >oo 



Therefore, when summing over k in evaluating the condensed series according to Eq. (|9|), it is in 
many cases sufficient to to evaluate the condensed series by adding the terms successively, and 
no further acceleration of the convergence is required. 

As shown in [|,15|, the condensation transformation defined according to Eqs. (p[)-(p!o|) is essen- 



tially a reordering of the terms of the input series X^^q a{k). Furthermore, Daniel was able to 
show (see the Corollary on p. 92 of Ref. fl^), that for nonalternating convergent series whose 
terms decrease in magnitude (|a(A;)[ > \a{k + 1)|), the equality (^) holds. This formally justifies 



the correctness of the condensation transformation defined according to Eqs. (|8|) - (10). 
Note that the property, originally derived in |p^], 

A2J-1 = \ (A,_i - a,_i) , (j = 1,2, . . . ) , (14) 

facilitates the numerical evaluation of a set of condensed series, by reducing the evaluation of 
condensed series of odd index to a trivial computation. Within the program, we use this relation 
in the form 

Aj+i = ^ (Ai/2 - ai/2) (even i) (15) 

in order to compute condensed series Am of odd index m = i + 1 where i is even. 

In the second step of the CNCT, the convergence of the Van Wijngaarden transformed series 
Z^jlo ("-'-)■' -Aj on the right-hand side of Eq. ^ is accelerated by a suitable nonlinear sequence 
transformation. We start from the partial sums 

n 

Sn = E(-l)'Ai (16) 

j=0 

of the Van Wijngaarden transformed series. Formal operations given by operators acting on the 
space of sequences are of considerable importance for the construction and theory of sequence 
transformations. We define the (forward) difference operator A 

A : {.4^=0 ^ {(As)„}^=o , (A5)n = Sn+l - Sn . (17) 

The nth element of the sequence {(As)„}J^o is given by (As)„ = s„+i — Sn- This way of writing 
the nth element of the transformed sequence stresses the fact that (As)„ is to be regarded as an 
element of some new sequence. However, the brackets are often left out, and we write 

As„ = Sn+l - Sn = (As)„ . (18) 

A well-known important relation is 

iA''s)n = i-l)'J2(-iy ^ j Sn+,. (19) 

As before we identify 



(20) 



The kth power of the difference operator plays a crucial rule in constructing sequence transfor- 
mations [see e.g. 0]]. A very important class of sequence transformations ("Levin-type transfor- 
mations"), which are characterized by the use of explicit remainder estimates is constructed 
according to the following "prescription" 0,|5|,16,17] 



cHfo N _ A'' {{n + (3)k-l Sn/uJn} 
<->k [PTSn,UJn) — — A J. r/ 7 ^ 





{P + n + j)k-l Sn+j 

j j {f3 + n + k)k-i iOn+j 

(21) 



k \ {(3 + n + j)k-i 1 



j j {(3 + n + k)k-^i iVn+j 

(n) 

The quantities appearing in the definition of 5^ have the following interpretation: 

• Sn- elements of a sequence whose convergence is to be accelerated. 

• ujn- remainder estimates, i.e. estimates for the truncation error ujn ~ rn defined in (^). 
Various sequence transformations based on Eq. (|2^), which differ only in the choice of 
uj-n, have been discussed in detail in [^. For all these sequence transformations, iOn can 
be calculated on the basis of the Sn- E.g., for the delta transformation defined below in 
Eq. ([2^), which is exclusively studied in the following, the calculation of tOn = Asn = 
Sn+i — Sn requires as input the sequence elements s„ and Sn+i- 

• l3: a numerical shift parameter, which may be adjusted or optimized for a particular appli- 
cation, e.g. by establishing asymptotic relations between the uJn and the r„, or by numerical 
experimentation. In practice, it is observed that any choice other than /? = 1 does not lead 
to an appreciable improvement of the rate of convergence of the transforms sj^^ , and there- 
fore f3 is almost exclusively set to unity in the literature (see also the numerical examples 
in [Q, ^). In order to allow for the case n = in Eq. (^), we must demand that /3 7^ 1, 
and in order to ensure a regular behaviour of the remainder estimates as a function of the 
index n + j, we should choose /3 > 0. 

• k: order the transformation. The forward difference operator A acts on the index n. The 
fcth power of that difference operator enters into the numerator and into the denominator 
of Eq. (pi). 



(n) 

• n: "initial element" of the transform. The computation of the transform Sf, requires 
as input the k + 1 elements . . . ,Sn+k} of the sequence {sn}^=Q, together with the 
k + 1 remainder estimates {cOn, ■ ■ ■ ,uJn-\-k}- The element of lowest index required for the 

(n) 

calculation of is Sn, and therefore n can be interpreted as the "initial element" in the 
calculation of the transform 5^"^ . 

At present, no mathematical proof exists with regard to the convergence of the transforms sj^\ 



although the construction of sequence transformations of the type (21) can be motivated heuris- 
tically (by expanding the truncation error in a factorial series) and is supported by overwhelming 
numerical evidence |p|,^,p^,p^. 

The notation S^\f3, Sn,uJn) is in need of a certain further explanation. As explained above, n 
represents the initial element of the sequence {s„}J^q used in the evaluation of the transform. 
The specification of Sn as the second argument of sj^^ (/5, Sn, Wn) is therefore redundant as far as 



the index n is concerned, because n already appears as an upper index in sj^"^ . The arguments 
Sn (and LVji) are to be interpreted as follows: they rather specify the initial element Sn and the 
initial remainder estimate ujn which are required for the calculation of the right-hand side of 
(|2l|). Of course, the evaluation of sj^\f3, Sm^Jn) requires in total the k + I sequence elements 
{sn, ■ ■ ■ , Sn+k} and the k + 1 remainder estimates . . . , uJn+k} as input. The knowledge of the 
"post-initial" sequence elements {sn+i, ■ ■ ■ , Sn+fc} and remainder estimates {un+i, • • • , ^n+k} is 
by assumption guaranteed in writing the definition (|2l|). A somewhat more clumsy notation for 



sj^\(3, Sn,uJn) which avoids all redundancy would read 

^k^\P^ {sjj^O' {^}j^o) (alternative notation not used in this article). 

The delta transformation @, ^] is a special case of the transformation (|2l|). It is constructed 
according to 

4")(/3,s„) = 5f)(/3,s„,A.O . (22) 



We use as input data for the delta transformation the partial sums (|16|) of the Van Wijngaarden 
transformed series. This corresponds to the "replacements" 

Sn^Sn, As„^AS„ = (-!)"+! A„+i, (23) 

and here, we always use So as the initial element for the second transformation (see however the 
discussion on p. 203 of where it is indicated that some terms of a given input series may 
have to be skipped in practical applications where the first terms of the input series display a 
highly irregular behaviour). This leads to the CNC transforms 

TcNc(n) = 4°)(l,So) , (24) 

which require as input the elements {So, . . . , Sn, Sn+i} of the Van Wijngaarden transformed 
series. 

A final word on sequence transformations: It is observed that sequence transformations con- 
structed according to ( ^l]) often lead to much better numerical results than better known con- 
vergence accelerators such as the Euler transformation, Wynn's epsilon and rho algorithms, 
or Aitken's process (all of these algorithms are described in [Q]). In view of the partially 
unsatisfactory situation concerning the availability of mathematical proof, we have carried out 
high-precision studies of the convergence of sequence transformations of the type (|2l|) , in part 
using multi-precision arithmetic. In the absence of a rigorous proof, one of the concerns which 
may be raised against the convergence of the transforms (^) is a plausible asymptotic nature of 
the sequence of transforms, resulting in ultimate divergence in higher order. In our numerical ex- 
periments, we found no indication of unfavorable asymptotic behaviour. The rate of convergence 
of the transforms remained constant, and we observed apparent convergence up to 430 decimal 
figures (see Sec. |^ below). In the absence of rigorous proof and the presence of considerable 
numerical evidence for the computational usefulness of sequence transformations of the type 



(21), it is conceivable that "mathematical experiments" may give some hints at the theoretical 



soundness of the concepts involved. 



3 Applications to Special Functions used in Statistics 

Several slowly convergent series of the type considered in this article define special functions that 
have important applications in statistics. These functions are Riemann's zeta, the generalized 



zeta, Jonquiere's function, and Lerch's transcendent (see |20| for definitions and some properties 



of these functions). Discrete distributions are related to these functions by associating probability 



mass functions (p.m.f.s) with the terms of the infinite series that define the functions. In turn, 
the normaUzation constants of the p.m.f.s are associated with the sums of the series. 

Consider a discrete distribution with the p.m.f. 

Fr[X = n]= cp{n; 9), (25) 

and with support of all nonnegative integers n = 0, 1, . . . , where is a (vector of) parameter(s) 
and c is a constant. The probability over the set of all outcomes is unity: 

oo 

Y^PT[X = n] = l, (26) 

n=0 

from which we obtain the normalization constant 

oo 

c-' = ^p{n;9) (27) 

n=0 

(in statistics, the random variable is usually denoted by the symbol x, whether continuous or 
discrete, but we will use the symbol n here, which is the more common notation for the running 
index in the theory of special functions). Due to its very property as a normalization constant, 
we expect the sum in Eq. ( p7|) to exist and to converge to a finite value. In the following, we 
consider several related discrete distributions. The Zipf distribution has the p.m.f. 

Pn = C^ (28) 



with support of all positive integers n = 1, 2, . . . . The normalization constant (|27| ) is then 

oo ^ 

= Cis) = E ' (29) 

n=l 

where s > 1. This is immediately recognized as the Dirichlet series for Riemann's zeta function 
[Eq. (1) on p. 32 of Ref. |2^]. Basic properties of the Zipf distribution can be found in [21|. 

The Zipf-Mandelbrot distribution is the generalization of the Zipf distribution that adds a 
constant v to the ranks n and has the p.m.f. 

Pn = c 7 — , (30) 

with the support of all nonnegative integers n = 0, 1, ... . The normalization constant (^) is 
then 

oo ^ 

c-^ = C{s,v) = Y- -, (31) 

n=0 



where v ^ 0, —1, .... The sum ( |3l|) is recognized as the generalized zeta function [Eq. (1) on 
p. 24 of Ref. [0]. 

Further, the p.m.f. 

Pn = c— (32) 

with the support of all positive integers n = 1,2,... defines the Good distribution |23|. The 
normalization constant (|2^) is then 

oo „ 

c-' = Fiz,s) = Y,'^, (33) 

n=l 



where \z\ < 1, which is recognized as Jonquiere's function [Eq. (14) on p. 30 of Ref. IQ; see 
also [^]. Kulasekera and Tonkyn derived moments of the Good distribution and proposed 
parameter estimators (for parameter estimation see also p6|). 

Now we discuss briefly several applications of the above distributions. While Zipf and Zipf- 
Mandelbrot distributions have been traditionally used in linguistics and information science, as 
well as in other disciplines |22, 2^, they also appear as canonical distributions in the general- 
ized statistical-mechanics ensemble theory. A variety of physical systems e.g. with long-range 
interactions or fractal boundary conditions are intractable within the classical Boltzmann-Gibbs 
statistical mechanics, while they can be successfully handled within the generalization proposed 
by Tsallis |^]. In particular, the variation of the generalized entropy with appropriate constraints 
yields power-law equilibrium distributions of the Zipf-Mandelbrot type |29|. The convolution of 



these distributions was shown to converge to the limiting Levy distribution [30|, according to 



the Levy-Gnedenko generalized central- limit theorem |31]. Furthermore, the interrelationship 
between microscopic statistical and macroscopic thermal equilibria for such power-law systems 
was shown to be a natural extension of the corresponding Boltzmann-Gibbs theory |3^ . 

The Good distribution has been used in linguistics to model word frequencies |3^, 34|. It was 
originally formulated by Good by introducing a "convergence factor" z" in Zipf's p.m.f. ([2^ ) 
in order to make the first moment of the distribution finite for 1 < s < 2. The p.m.f. (p2|) as an 
approximation for the upper tail of the Yule distribution [^] was mentioned by Simon [p^ . The 
particular case of the Good distribution with s = 1, known as Fisher's logarithmic law |^^, has 
been obtained as an equilibrium solution of certain stochastic processes by Champernowne in 
his model of the distribution of incomes |37|, and by Darwin in his model of the frequencies of 
species in evolving populations |38|. A logarithmic law for the frequencies of species has also been 
obtained as a limiting form of the negative binomial distribution The Good distribution 
was also suggested for use in ecology for the reasons of the flexibility of its hazard function . 

In biophysics, one is interested in establishing the statistical structure of various sequences, e.g. 
those of DNA, RNA, or proteins. Indeed, an interest in statistical compositions of DNA and 
protein sequences originated shortly after the discovery that the nucleotide sequences of RNA 
uniquely determines the amino-acid sequence of the proteins. An early analysis was performed by 
Gamow and Yeas , who observed that the distributions of relative abundances of amino-acids 
and nucleotides were nonrandom, i.e. deviated from models assuming a uniform distribution of 
different words in any given text. Further studies established the utility of relative abundances 
of words (short oligonucleotides) in a comparative analysis of genome sequences [41,42|. Further- 
more, a linguistically motivated "Zipf" analysis has been applied to a number of DNA sequences 
in an attempt to demonstrate that noncoding sequences resemble a "natural" language more 
than the coding ones [43|. However, it was noted that the differences between the fitted Zipf 
exponents s in the p.m.f. (|2^) are not indicative of the closeness to a "natural" language per 



se 1 44 1, as "Zipf-like" (i.e., power-law) frequency distributions arise in a wide variety of situa- 
tions. Moreover, the obtained fits do not apply in the whole range of word ranks (i.e., slopes of 
empirical plots in log- log coordinates are not constant). In contrast, the Good distribution was 
found to fit the data in a wider range of ranks 1^^, even though it was not possible to determine 
the statistically significant differences between the fits to coding and noncoding sequences, due 
to large sampling errors of the fitted parameters. These studies reinforce an early comment by 
Simon |33] that the possible cause for the wide applicability of the "power-law" Zipf-related 
distributions is perhaps the similarity in the probability mechanisms underlying such diverse 
phenomena as distributions of words in DNA sequences and in "natural" languages, frequen- 
cies of citations of scientific articles and of species abundances in populations. Interestingly, the 
possibility of such a unification might occur within the framework of the generalized thermo- 
statistics, which is exemplified by obtaining the Zipf-Mandelbrot law in binary sequences within 
the generalized thermodynamics of Tsallis . 

A formal generalization of the Zipf, the Zipf-Mandelbrot, and the Good distributions is possible 



Distribution Support p.m.f. normalization 

(related function) constant 



Zipf 1,2,... cn-" Cis) = ^(1,8,1) 
(Riemann's zeta) 

Zipf-Mandelbrot 0,1,... c(n + u)-^ C{s,v) = ^{l, s,v) 
(generalized zeta) 

Good 1,2,... cz"n-^ F{z, s) = z <i>{z, s,l) 
(Jonquiere's function) 

Lerch 0,1,... cz"'(n + t;)"* ^{z,s,v) 
(Lerch's transcendent) 



Table 1: Relationships between distributions defined by Riemann's zeta, generalized 
zeta, Jonquiere's and Lerch's functions. 



upon realizing that Riemann's zeta, the generalized zeta, and Jonquiere's functions constitute 
special cases of Lerch's transcendent which is defined by the following series 

OO Yi 

$(z,s,t;) = V-^— , (34) 

n=0 ^ ' 



where \z\ < 1 and v ^ 0,-1,... [Eq. (1.11.1) on p. 27 of Ref. pO[]- '^^^ relations between 
$(z, s, v) and C{s), C{s, v), and F{z, s) can be obtained by letting z and/or v to and/or 1, and 
making use of the following functional relation [Eq. (1.11.2) on p. 27 of Ref. |20| ]] 



<^{z,s,v) = z"'^{z,s,m + v) + y2-——. (35) 

In particular, Eq. ( ^5[ ) with m = 1 can be used to rewrite sums that start from n = 1 rather 
than from n = in terms of the ^ function: 

OO „ 

Y.^^-— = z^z,s,l + v) (36) 

n=l 

In Tab. |l], we have summarized the relations between the different statistical distributions and 
their defining special functions. The three distributions, Zipf (zeta), Zipf-Mandelbrot (general- 
ized zeta), and Good (Jonquiere), are special cases of the Lerch distribution, and their properties 
can be expressed in terms of Lerch's transcendent with special values of the parameters. (See 
for the discussion of relations between these distributions and the zero-truncated Lerch distri- 
bution.) 

Here we consider several functions of the related Lerch distributions, from which other properties 
(e.g., moments) can be easily derived (see |47| for details on the Lerch distribution). 



The cumulative distribution function (c.d.f.) F{n; 9) of the Lerch distribution with the support 
n = 0, 1, . . . is given by the following equation 

F(n; 9) = Pr[X < n] = 1 - + ^ + . (37) 

The c.d.f.s for the Zipf, the Zipf-Mandelbrot, and the Good distributions are given by 

_ ^(l,g,re+ 1) _ $(l,g,^ + n+ 1) _ z"$(z,g,n+l) 

$(1,.,1) ' $(1,5,^;) ' ' ^""^^ 

respectively. The survival and hazard functions are obtained from the c.d.f.s ( |37| ) and ( p^ ) as 
5(n; ^) = 1 — F(n; 0) and /i(n; 0) = p{n; 9) /{I — F{n; 9)), respectively. 

The probability generating function (p.g.f.) G{y;9) for the Lerch distribution is given by 

G(y;^) = f;,"Pr[X = n] = ^^, (39) 

where we assume \y\ < 1. The p.g.f.s for the Zipf, Zipf-Mandelbrot and Good distributions are 
then given by 

y^iy,s,l) '^{y,s,v) y<^{yz,s,l) 

respectively. The characteristic functions (c.f.'s) and moment generating functions (m.g.f.s) can 
be obtained from the p.g.f.s ( |39| ) and ([40|) by letting y = e** and y = e*, respectively. The 
moments of the Lerch distributions are then obtained as coefficients of the Taylor series of the 
m.g.f.s about t = 0. 

We note that so far the support of the distributions consisted of all nonnegative integers, except 
possibly 0, whereas in applications the distributions are often used in singly or doubly truncated 
forms. Truncation of a distribution changes the way the normalizing constant is calculated. 
For the truncated case of the Lerch distribution with support n € [a,b], where a > and 
a <b < 00, the c.d.f. is now 



n 



k 



z°- <l>(z, s, w + a) - 2;"+^ $(z, s, f + n + 1) 



^ ' ^ ^(k^vY z'^4>(z,s,t; + a) -z''+i^>(z,s,?; + 6+l) ■ ^ ^ 

k=a 

Likewise, the p.g.f. for the Lerch distribution is 

Q( .0)- (y ^iyz,s,v + a)- {y z)^+^ ^{yz, s,v + b+l) 

For truncated Zipf, Zipf-Mandelbrot, and Good distributions, the c.d.f.s can be obtained from 
Eq. (PI), and the p.g.f.s can be obtained from Eq. (p2|), by substituting z = 1 and v = (for Zipf), 
z = 1 (for Zipf-Mandelbrot) and v = (for Good). Note that these "short-cut" substitutions 
are possible because truncated distributions have the same support (by definition). We cannot 
use these substitutions to obtain Eq. (|^) and (|40|) because in that case the distributions have 
different support, depending on whether zero is excluded or included (see Tab. ||). 

Calculations with the Lerch distributions using Eqs. (^) - (42) involve an evaluation of Lerch's 
transcendent. Other commonly needed manipulations with the Lerch distributions, e.g. random 
number generation and parameter estimation, can be likewise reduced to calculating the Lerch 
transcendent [HtII. 



We have implemented the calculation of Lerch's transcendent ^{z,s,v) defined in Eq. (34) 
using the CNC transformation [Eqs. (|8|) - (p!oD] in Mathematica. The program is presented 
and discussed in Appendix For a provisional implementation in C as well as a Mathematica 
package for calculations with the Lerch distribution, we refer to 



4 An Application in Experimental Mathematics 



We begin this section by quoting [10|: "In April 1993, Enrico Au-Yeung, an undergraduate at 
the University of Waterloo, brought to the attention of the author's [David Bailey's] colleague 
Jonathan Borwein the curious fact that 



E 

k=l 



1 + 5 + 



4.59987- 



17^ 
360 ' 



(43) 



based on a computation of 500,000 terms. Borwein's reaction was to compute the value of this 
constant to a higher level of precision in order to dispel this conjecture. Surprisingly, his compu- 
tation to 30 digits affirmed it. [David Bailey] then computed this constant to 100 decimal digits, 
and the above equality was still affirmed." Many formulas similar to have subsequently 



been established by rigorous proof [5C]. 

With the help of a multiprecision system [|^10j and the CNCT, we have verified ( |43| ) "experi- 
mentally" to a couple of hundred decimals [51]. The calculation will be sketched in the following. 
Using the definition 



1+-+...+.I k- 



(44) 



we rewrite 



as follows, 



OO OO / fc \ 

k=Q k=0 \j=0 I 



E 

A:=0 



fc-h 1 



(45) 



five of the Gamma function 1 2C , 5S 



where 7 is the Euler-Mascheroni constant 7 = 0.577 215 . . . , and i/j^z) is the logarithmic deriva- 

(46) 



With the help of the relation 



fc=o 



E|^-c(3)-.c(.,, 



Eq. (Hsl) can be rewritten as 



E 

fc=0 



V'(fc + 2 )V _ 17 
k + 1 



■C(4)-47C(3)+7'C(2). 



We proceed to calculate numerically, to high precision, the infinite sum 



(47) 



(48) 



(49) 



using the CNC transformation. 



In order to establish the rate of convergence of (^), we investigate the asymptotic behaviour of 
the a{k) as k ^ oo. The logarithm of the Gamma function can be expanded into an asymptotic 
series [see Eq. (4.03) on p. 294 of 



Inr(z) = (^z--y^z - z + - ln(27r) + ^ 2.(2. -l)z^s-i + ^-(^) ' (^0) 



where 



Here, \x\ is the integral part of x, i.e., the largest integer m satisfying m < x, Bk[x) is a 
Bernoulli polynomial defined by the generating function [see Eq. (1.06) on p. 281 of Ref. |5^ ]: 



B2m{x - {x 



m=0 

and 



tex.r>(xt) V— ^ t"^ 

i5;(|zT = I'l<2-. (52) 



Bm = B^{0) (53) 
is a Bernoulli number (p. 281 of Ref. Is^). The following asymptotic relation for ^^(-2) follows: 



S = l 



Z^ 



The leading asymptotics of the remainder of the sum ( ^9|) after adding — 1 terms can thus be 
derived easily. We have for large k, 

ln(fc + 2f ln{k + 2) ln(fc + 2) , f I \ , 

"(^^ ^ 7feTTF" (fc + l)Mfe + 2) " 6(fc + l)Mfc + 2)^ +^VW ' ^''^ 



Based on these formulas, the remainder of the sum (|49|), for large N, can be written as 

k=N ^ ^ 



Here, the Euler-Maclaurin formula [Eqs. (2.01) and (2.02) on p. 285 of Ref. [52|] has been used 
in order to convert the sum over the a{}i) in the asymptotic regime of large k [see Eq. (|55|)1 
into an integral plus correction terms. In order to calculate (^) to an accuracy of 200 decimals, 
Eq. ( |56[ ) says that we would be required to add on the order of 10^'^^ terms. Without the use of 
convergence acceleration methods, this would represent a formidable computational task. 



Using the CNCT, it is easy to calculate the sum ( |4q ) to 200 digits, based on multiprecision 
arithmetic [^10| and a Linux personal computer, within a few hours. We obtain for the 246th 
and the 247th CNC transform defined according to Eq. (24), 

rcNc(246) = 2.37254 51620 38445 67035 68130 69148 85258 25756 18499 54254 



97013 57806 20011 72404 62937 46020 32218 23862 67095 00004 



69194 36541 28946 10390 15116 52595 90270 23975 58737 74256 



23420 48480 95165 00802 19816 35378 76591 98589 60393 32102 8 , (57) 



and 



TcNc(247) = 2.37254 51620 38445 67035 68130 69148 85258 25756 18499 54254 
97013 57806 20011 72404 62937 46020 32218 23862 67095 00004 
69194 36541 28946 10390 15116 52595 90270 23975 58737 74256 

23420 48480 95165 00802 19816 35378 76591 98589 60393 32111 7. (58) 

The apparent convergence to roughly 200 decimals can be verified against the right-hand side 
of Eq. (^8|). Of course, the right-hand side of Eq. (|48|), which involves only rationals, zeta 
functions and the Euler-Mascheroni constant, 

^C(4)-47C(3) + 7'C(2), 

can be easily evaluated to 200 decimals using known algorithms which are included in computer 
algebra systems (e.g. |^,^)- 

The evaluation of the terms a(fc) proceeds as follows. For small index it is easy to write a 
recursion relation relating a{k) and a{k + 1) based on the (trivial) recursion for the ip function, 

il>{k + l)=^{k)+'^. (59) 

For large A;, the asymptotic formula (|5^ ) can be used in order to calculate the ip function to 
high precision. The point at which one may switch from the recursion to the asymptotic method 
depends on how many explicit values for Bernoulli numbers are available to the machine. We 
use the values for the first 60 Bernoulli numbers, to 250 decimals, for our calculation. We switch 
from one method to the other when the index k of a{k) has reached a value of 500. 

With 84 308 a{k) terms evaluated (out of which 1364 by recursion and 82944 by the asymptotic 
method), we evaluate the first 247 transforms with the results presented above in Eqs. (|57| ) and 
(|58D. As mentioned above, if the terms of the series ( |4^ ) were added on a term-by-term basis, then 
about lO^''^ would be required for an accuracy of 200 decimals in the final result. The reduction 
of this number to roughly 84 000 corresponds to an acceleration of the convergence by roughly 
200 orders of magnitude. We have also carried out, using enhanced precision, a calculation to 430 
decimals, involving about 500 CNC transformations and arithmetic with 600 decimal figures. 
These evaluations not only confirm the relation (|48| ) to high precision, but they also represent an 
accurate experimental verification of the convergence properties of the delta transformation ([2^ ) 
in higher transformation order. Moreover, it is observed that the rate of convergence of the CNC 
transform results in a gain of approximately one significant decimal figure per transformation, 
and that the rate of convergence remains constant over a wide range of transformation orders 
(corresponding to linear convergence in the asymptotic region). By contrast, the series ( ^8|) is 
only logarithmically convergent. This corresponds to convergence acceleration according to the 
definition (^). 

5 Other Applications of the CNCT 

We also mention the existing applications in the domain of quantum electrodynamic bound-state 
calculations (see |2|,|l^,^. Another application concerns the quantum electrodynamic effective 



action (see [^). The combined nonlinear-condensation transformation is also applicable |51| to 
series of the form 

oo 2fe+l 

R jx) = y -4 — ^, (60) 



( h\ = 1 cosh(2fc + l)x 

^ ^ (2A; + 1)P cosh(2fc + 1)6 ' ^ ' 



fc=0 



and 

Up{x 



M = V 1 cosh(2fc + 1)3; 

^ (2A; + 1)P sinh(2A: + 1)6 ' ^ ^ 



Series of this type occur naturally in the context of plate contact problems with periodic bound- 
ary conditions ||5^. The arguments p, x and 6 are real and positive for cases of practical rele- 
vance. For X ~ 6 and p ~ 1, the series Tp and Up are very slowly convergent. In App. A. 2 of 
(pp. 105 ff. ibid.), it is demonstrated by way of numerical example that the CNCT is able to 
efficiently accelerate the convergence of these series in problematic parameter regions. 

In the numerical calculations, it is necessary to evaluate terms with large index k. This can 
lead to numerical overflow because of the large arguments of the hyperbolic functions. Clearly, 
representations such as 

Tp{x,b) = (2fc + l)P 1 + e-2K2fc+i) 

fc=0 

provide a solution for this problem. 

Let us recall that considerable effort has been invested in the development of efficient numerical 
methods for the evaluation of the series ( |60| ) - ( |62[ ) |55-59|. These alternative methods make 
intensive use of special properties of the series. They involve integral transformations and infinite 
series over numerical integral [^, and they make use of special properties of Legendre's chi- 



Function which is related to the functions (pOD - (|6: 



We also briefly mention that it is possible, in combining analytic results obtained in [11, 60|, 61 
with numerical techniques based on the CNCT for the evaluation of complete hypergeometric 
functions, to evaluate the so-called Bethe logarithm in hydrogen to essentially arbitrary precision. 
Specifically, we obtain - for the 4P state - the result 

lnfeo(4P) = -0.041 954 894 598 085 548 671 037(1) , (64) 

which should be compared to other recent calculations ||6^|6^. 



6 Conclusions 

We have discussed several applications of the convergence acceleration methods introduced in 
Sec. |2|: in statistical physics (Sec. |3|), in experimental mathematics (Sec. and other appli- 
cations, mainly in the evaluation of special functions (Sec. ^). Specifically, it is observed that 
the combined nonlinear-condensation transformation (CNCT) leads to an efficient calculational 
scheme for the Lerch transcendent $ given by Eq. (|3^. This special function provides a general- 
ization of several kinds of probability density functions which are of significance for the statistical 
analysis e.g. of DNA sequences (see Table ^). The comparatively fast and accurate evaluation 



of the statistical distributions with the CNCT is helpful with respect to statistical procedures 
such as random number generation and statistical inference |47|. 

The high-precision calculation described in Sec. ^ confirms that the rate of convergence of the 
CNC transforms is consistent with linear convergence (see Sec. |2.1| ) in large transformation 
orders, and in extended numerical precision (up to 430 digits). In the absence of a rigorous 
proof Q regarding the convergence of the nonlinear sequence transformation (22), which forms 
the second step of the CNCT (see Sec. 2.2), we attempt to gain numerical evidence for the 
(linear) convergence of the CNCT in higher transformation orders, and, in Sec. ^, for the wide 
applicability of this algorithm in the context of the special functions which are relevant to 
theoretical physics. 
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A Mathematica Program for Lerch's Transcendent 

The listing of the Mathematica program for the calculation of Lerch's transcendent is shown 
in Fig. |l|. The code displayed in Fig. || is compact, but it lacks input argument checking. Also, 
there is no convenient control over the following parameters: (i) relative accuracy of the result, 
(ii) maximum number of allowed iterations, (iii) shift parameter P, and (iv) "initial element" 
parameter n. The code should be understood as a template for more specialized implementations. 
The program in Fig. |l] takes advantage of Mathematica's implementation of the functional 
programming paradigm. 

Necessarily, the Mathematica code shown in Fig. |l] is restricted to positive argument z > 0. 
For the other arguments, the same restrictions exist as for the C program available from [pS[|: 
s,v must be real, the case of negative integer v is excluded, and for negative non-integer v, s is 
required to be integer. 

The overall structure of the program is the following: lines 2-19 contain declarations of the 
three subroutines that are used iteratively during the calculation, and lines 20 - 24 contain the 
main body. We explain the algorithm underlying the program on a line-by-line basis. 



1: Declare the LerchPhiCNCT function. 

2: Declare the subroutine ajstep that performs one iteration for calculation of the Van Wi- 
jngaarden transforms according to Eqs. (|^ and (10). The subroutine ajstep comprises 
the program lines 2-5. 

3: Read the current iteration counter k and the temporary value of the partial sum of the 
series defining Aj [see Eq. (^)] into local variables. 

4: Calculate a new , update the temporary value for the sum Aj, and increment the 
counter k. 



5: Return the result of one ajstep iterate, which is a list of the counter k, the current value 

(i) 

of partial sum defining Aj, and the last calculated . 



1 : LerchPhiCNCT [z_ , s_ , v_] : =Module [{j , omegaO, storeaj , nviin, den, o, res .factors} , 
2: ajstep[inp_List] :=Module[{k,sum,bjk}, 
3 : k=inp [ [1] ] ; siain=iiip [ [2] ] ; 

4: bjk=2-k z-(2-k (j+l)-l)/(v+2-k (j+l)-l)-s;sum+=bjk;k++; 

5: Return [{k, sum, b j k}] ] ; 

6: recur [inp_List ,pos_] :=Module[{loc}, 

7: loc=inp; loc [ [pos] ] =loc [ [pos+1] ] -loc [ [pos] ] f actors [ [pos] ] ; 
8 : Return [loc] ] ; 

9: sknstep [inp_List] : =Module [{i , omega, sn, skn,eps}, 

10 : i=inp [ [1] ] ; omega=inp [ [2] ] ; sn=inp [ [3] ] ; skn=inp [ [4] ] ; eps=inp [ [5] ] ; 
11 : sn+=omega; AppendTo [storeaj , (-1) "i omega] ; 

12: omega=(-l)~(i+l) If [EvenQ [i] ,0 . 5 (storeaj [[i/2+l]]-z-(i/2)/(v+i/2)-s) , 

j=i+l ;NestWhile [ajstep,{l ,z''j/ (v+j) ~s,z"j/(v+j) "s}, 

(#[[3]]/#[[2]]>10-(-2-acc))&] [[2]]] ; 
13 : AppendTo [num , sn/ omega] ; AppendTo [den , 1/ omega] ; 
14 : f actors=Reverse [Which [i==0 , {0} , i==l , {0 , 1} , True , 

Prepend [Table [(beta+n+i-1) (beta+n+i-2)/(beta+n+i+o-2)/(beta+n+i+o-3) , 

{o,l,i>] ,0]]] ; 

15: num=Fold[recur, num, Table [o , {o ,Length[num] -1 , 1 , -1}] ] ; 
16: den=Fold [recur, den, Table [o,{o, Length [den] -1,1,-1}]] ; 

17: skn=RotateLef t [skn] ; eps=RotateLef t [eps] ; 

18 : skn [ [2] ] =num [ [1] ] /den [ [1] ] ; eps [ [2] ] =Abs [skn [ [2] ] -skn [ [1] ] ] ; i++ ; 
19 : Return [{i , omega , sn , skn , eps}] ] ; 
20: acc=14;imax=100;beta=l;n=0; 

21 : j=0 ; omegaO=NestWhile [aj step,{l ,z" j/(v+j ) "s ,z~ j/ (v+j ) ~s} , 

(# [ [3] ] /# [ [2] ] >10~ (-2-acc) ) &] ; num={} ; den={} ; storeaj ={} ; 
22 : res=NestWhile [sknstep , {0 , omegaO [ [2] ] , , {0 , 0} , {0 , 0}} , 

Not[#[[l]]>imax-l| I (#[[1]]>1&&(#[[5,2]]==0| I (#[[5,2]]<#[[5,1]]&& 

Abs[2 (#[[5,l]])-2/(#[[5,l]]-#[[5,2]])/#[[4,2]]]<10-(-acc))))]&] ; 
23: If [res[[l]]>imax-l. 

Print ["Algorithm has not achieved relative accuracy of ",acc, 

" digits after maximum "res[[l]]," iterations."]]; 
24: Return [res [[4, 2]]]] 

Figure 1: A Mathematica program for the calculation of Lerch's transcendent using the CNC 
transformation. A line- by-line explanation is in the text. 



6: Declare the subroutine recur that performs one iterative recursive evaluation of the numer- 
ator and the denominator of the sequence transformation in Eq. (^T]). The algorithm takes 
advantage of the recursion relation in Eq. (3.11) of Q. The subroutine recur comprises 
lines 6 - 8 of the program. 

7: Read the input into a local variable and apply the recurrence relation. 

8: Return the result of recur, which is an updated list loc. 

9: Declare the subroutine sknstep which performs iteratively the calculation of the CNC 
transforms. The subroutine recur comprises lines 9 - 19. 

10: Read the following data into local variables: (i) the current iteration counter i, (ii) the Van 
Wijngaarden term (also remainder estimate) uji-i = (— 1)* Aj, (iii) the partial sum of the 
Van Wijngaarden series Si (which is the input data for the delta transforms), (iv) the last 
two CNC transforms calculated [6l^2 = ™^ '^l-il' i^) difference between CNC 
transforms of successive iterations ej_2 and ej_i where = \6^^^ — The definitions 

of these symbols can be found in Eqs. (|) - (|^), (||), (||), and (H). 

11: Increment the partial sum Sj with the next term (—1)* Aj of the Van Wijngaarden series, 
and store the current Aj in the list of variables storeaj for later use in the recurrence 
relation Eq. ([l5|). 



12: Calculate the next term uji = (—1)*"^^ Aj+i of the Van Wijngaarden series: for even index 
i, use Eq. ([I5|). For odd i, calculate via a call to the subroutine ajstep. The last argument 
in the "nested call" to the Mathematica built-in function NestWhile is a termination 
criterion: the are added on a term- by-term basis until the ratio of the third element 

of the list of the results of ajstep (which is the value of the current b^"'^) and the second 
element of that list (which is the value of the current partial sum of Aj) is less than 
the specified accuracy. This termination criterion is justified for the case of the Lerch 
transcendent by Eq. (13). 



13: Calculate and append new values in the lists that store the numerator and denominator 
of the 6f^^sf\ 
transformation i. 



of the df'^ = sji^^ . The length of these lists corresponds to the current order of the CNC 



14: Calculate a list of prefactors used for the recursive evaluation of the denominator and 
numerator of the delta transformations df'^ according to Eq. (3.11) of 

15: Calculate the numerator of the delta transforms recursively and store the results in the 
list num. 

16: Calculate the denominator recursively and store results in the list den. 

17: Swap the elements of the lists skn and eps which contain the last two delta transforms 
and the last two absolute differences; the purpose is to let the "old" values be the first 
elements, in order to store the "newly calculated" values as the second elements in a later 
step. 

18: Calculate a "new" delta transform and a new absolute difference eps and store them as 
second elements of skn and eps. 

19: Return the result of sknstep. It consists of a "new" iteration counter i, a "new" Van 
Wijngaarden term omega, a "new" partial sum of the Van Wijngaarden series sn, a list 
skn of the last two CNC transforms, and a list eps of the last two differences between 
CNC transforms of successive iterations. 



20: Start the main program. Set the default values for the desired relative accuracy of the 
result, for the maximum allowed order of the delta transformation, for the shift parameter 
/3 = 1, and for the "initial element" parameter n = which enters the CNC transform 
in Eq. (H). 

21: Calculate the initial remainder estimate uj via a nested call to a j step (in analogy to line 
12) and initialize other variables. 

22: Iteratively compute the CNC transforms using sknstep. Again, in analogy to line 12, 
the second argument of the NestWhile function is the list of initial values, and the third 
argument is the termination criterion. The evaluation of the CNC transforms stops if either 
the maximum number of iterations is exceeded (first element of the output of sknstep), 
or if the current difference between the transforms is zero (fifth element), or if an error 
estimate for the last evaluated delta transform is less than the specified accuracy. This error 
estimate is inspired by the fact that the delta transforms are observed to converge linearly 
in good approximation [see Eq. (|^ and Sec. |2.1| as well as Sec. |^. The remainder estimate 
is evaluated using the ej_i and (calculated from the fifth argument to sknstep) and the 
current delta transform (which is second element of the fourth argument to sknstep). 

23: Print an error message if the maximum number of iterations is exceeded. 

24: Return the result of LerchPhiCNCT, which is the second element of the fourth argument 
of sknstep and corresponds to the last delta transform evaluated. 



To conclude this Appendix, we provide a more detailed explanation for the error estimate used 
in program line 22. 

The output of the CNC transformation is a sequence of approximants Tn = '7cNc('^) [see Eq. (|2^)] 
that converge to the value of ^{z,s,v). We define the ratio of two consecutive differences of 
approximants as 



T - T 



J. ri 



n — 1 



(65) 



As discussed in Sec. ^ the CNC transforms can be expected to converge geometrically in higher 
transformation orders. Let p be defined as [see also Eq. ^\ 



lim -5— 



(66) 



with Tcx> = ^{z,s,v). Clearly, the Xn will provide a good estimate for p at large n. A good 
estimate for the truncation error 7^ — T^o can thus be obtained by summing the geometric series 
Ylk^i I'^n — %i-i\ where our best estimate for pispK, x„. We may therefore use the following 
convergence criterion to terminate the calculation of the CNC transforms: 



2 



1 



T — T 



n-1 



Tn 



< acc 



(67) 



Here, acc is the specified desired relative accuracy of the result. The factor 2/xn in ( |67[ ) is a 
heuristic "safeguard factor" introduced with the notion of avoiding a premature termination of 
the calculation of successive transforms in the problematic case of two consecutive transforms 
accidentally assuming values very close to each other. Such a situation may arise before the 
asymptotic, geometric convergence sets in. The term in square brackets in (|67| ) represents the 
remainder estimate based on the geometric model Xn ^ p [see Eq. ([65|)]. Rewriting (|67|), we 
obtain the termination criterion used in line 22 of the program in Fig. ffl. 



Finally, we would like to point out that the above implementation relies on only one algorithm 
(CNCT), and that one cannot expect to obtain optimal performance in all parameter regions, 
let alone analytic continuations for those cases where the power series (34) diverges. As regards 
the evaluation of special functions with very large (excessive) parameter values, it is known 
that asymptotic expansions can provide optimal methods of evaluation (see, e.g. |^). These 
are not implemented in the code in Fig. |l|. For negative 2, in our case some improvement can 
be obtained by accelerating the rate of convergence of the resulting alternating series by the 
direct application |^, 49| of delta transforms to the power series (|3^; this avoids the potentially 
time-consuming Van Wijngaarden step (P) - (10) of the CNCT. However, the more problematic 
case of positive z, which results in a nonalternating series, requires the (full) combined nonlinear- 
condensation transformation implemented by the algorithm in Fig. |^. 
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